Interpretabilidad en modelos de aprendizaje automático

La interpretabilidad se refiere a la capacidad de traducir conceptos abstractos del modelo a una forma comprensible: ¿cómo se generan las predicciones?.

• El objetivo es extraer conocimiento relevante del modelo sobre las relaciones en los datos o aprendidas por el modelo.

• Modelos más interpretables permiten al analista entender mejor las predicciones, lo que es clave en economía y negocios (p.e., en finanzas, para la concesión de un préstamo, es relevante entender cómo el modelo genera puntuaciones crediticias con el propósito de evaluar la solvencia del cliente).

La interpretabilidad es relevante ya que nos permite identificar asociaciones entre predictores y la respuesta (que nos pueden llevar a formular relaciones causales) así como detectar y corregir problemas (p.e. sesgos) en el modelo estimado.

Fuente: An introduction to statistical learning

Fuente: Interpretable Machine Learning (Molnar, C.)

Modelos paramétricos interpretables

# carga de bibliotecas y datos 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

url = "https://juandmontoro.github.io/bigDataEco/data/housing.csv"
data = pd.read_csv(url)

# Seleccionamos los nombres de las var numéricas
var_numericas = data.drop(columns='median_house_value').select_dtypes('float64').columns
# Convertimos la var categórica a dummy
data = pd.get_dummies(data,dtype=int)
# Definimos matriz de diseño X y respuesta
X = data.drop(columns=['median_house_value'])
y = data['median_house_value']/1000  # expresamos el precio en 1000 USD
# Partición de datos
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=42)

# Instacia del estimador regresión lineal y ajuste
lm = LinearRegression()
lm.fit(X_train,y_train);

# Visualizamos los coeficientes
pd.DataFrame({'Variable':lm.feature_names_in_,'Coeficiente':lm.coef_})

                      Variable  Coeficiente
0                    longitude   -42.599568
1                     latitude   -41.977692
2           housing_median_age     0.360641
3                  total_rooms     0.050232
4               total_bedrooms    -0.176841
5                   population    -0.061597
6                   households     0.109672
7                median_income     0.003973
8    ocean_proximity_<1H OCEAN    -1.290712
9       ocean_proximity_INLAND   -71.612767
10      ocean_proximity_ISLAND    95.530903
11    ocean_proximity_NEAR BAY    -8.183977
12  ocean_proximity_NEAR OCEAN   -14.443448

from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer

# Preprocesamiento: incluimos transformaciones para var numéricas
ctr = make_column_transformer((make_pipeline(PolynomialFeatures(2,include_bias=False),StandardScaler()),
                               var_numericas),
                               remainder='passthrough')
# Estimador compuesto: preprocesamiento + estimador (regresión lineal)
pipe = make_pipeline(ctr,lm)
# Ajuste
pipe.fit(X_train,y_train);

# Visualización de coeficientes
# Necesitamos extraer los nombres de las variables creadas. Seleccionamos los elementos de pipe que transforman las variables
vars = pipe[0].get_feature_names_out()
pd.DataFrame({'Variable':vars,'Coeficiente (estandarizado)':pipe[1].coef_})

                                       Variable  Coeficiente (estandarizado)
0                           pipeline__longitude                  1178.500046
1                            pipeline__latitude                  2565.371678
2                  pipeline__housing_median_age                 -1010.716726
3                         pipeline__total_rooms                  -867.094149
4                      pipeline__total_bedrooms                  1825.479208
5                          pipeline__population                  2215.547121
6                          pipeline__households                 -3101.972606
7                       pipeline__median_income                   273.330277
8                         pipeline__longitude^2                  2214.274787
9                  pipeline__longitude latitude                  4466.065572
10       pipeline__longitude housing_median_age                 -1441.699979
11              pipeline__longitude total_rooms                 -1509.704720
12           pipeline__longitude total_bedrooms                  2321.035242
13               pipeline__longitude population                  3422.546382
14               pipeline__longitude households                 -4118.087702
15            pipeline__longitude median_income                   416.555707
16                         pipeline__latitude^2                   926.957257
17        pipeline__latitude housing_median_age                  -469.057414
18               pipeline__latitude total_rooms                  -519.579739
19            pipeline__latitude total_bedrooms                   411.366239
20                pipeline__latitude population                  1173.406307
21                pipeline__latitude households                 -1040.604400
22             pipeline__latitude median_income                   115.117677
23               pipeline__housing_median_age^2                    29.552256
24     pipeline__housing_median_age total_rooms                    42.716036
25  pipeline__housing_median_age total_bedrooms                   -14.285175
26      pipeline__housing_median_age population                   -47.964444
27      pipeline__housing_median_age households                    39.949436
28   pipeline__housing_median_age median_income                     2.624768
29                      pipeline__total_rooms^2                   336.535103
30         pipeline__total_rooms total_bedrooms                  -661.472108
31             pipeline__total_rooms population                  -396.387657
32             pipeline__total_rooms households                   309.294231
33          pipeline__total_rooms median_income                   -40.418864
34                   pipeline__total_bedrooms^2                   211.262862
35          pipeline__total_bedrooms population                   453.752989
36          pipeline__total_bedrooms households                  -107.943168
37       pipeline__total_bedrooms median_income                    60.426682
38                       pipeline__population^2                    67.060717
39              pipeline__population households                  -168.663219
40           pipeline__population median_income                    10.028064
41                       pipeline__households^2                   -18.435552
42           pipeline__households median_income                   -32.173088
43                    pipeline__median_income^2                    34.806495
44         remainder__ocean_proximity_<1H OCEAN                     7.968695
45            remainder__ocean_proximity_INLAND                   -63.818280
46            remainder__ocean_proximity_ISLAND                    67.726284
47          remainder__ocean_proximity_NEAR BAY                    -6.067859
48        remainder__ocean_proximity_NEAR OCEAN                    -5.808840

from sklearn.linear_model import LassoCV
lambdas = np.logspace(-3,4,200)
pipe = make_pipeline(ctr,LassoCV(alphas=lambdas,cv=5,max_iter=10000))
pipe.fit(X_train,y_train);
coefLasso =pd.DataFrame({'Variable':vars,'Coeficiente (estandarizado)':pipe[1].coef_})
# Filtramos el DF para eliminar las filas con coeficientes nulos
coefLasso[coefLasso['Coeficiente (estandarizado)']!=0]

                                    Variable  Coeficiente (estandarizado)
0                        pipeline__longitude                   -74.647785
5                       pipeline__population                   -34.544223
11           pipeline__longitude total_rooms                   -88.415017
15         pipeline__longitude median_income                    13.766031
16                      pipeline__latitude^2                   -70.570682
17     pipeline__latitude housing_median_age                   -48.330092
18            pipeline__latitude total_rooms                     4.379371
19         pipeline__latitude total_bedrooms                   -70.648855
21             pipeline__latitude households                    -6.102056
22          pipeline__latitude median_income                    -1.067559
23            pipeline__housing_median_age^2                    35.689047
24  pipeline__housing_median_age total_rooms                    42.873616
26   pipeline__housing_median_age population                   -36.598430
27   pipeline__housing_median_age households                    19.948919
30      pipeline__total_rooms total_bedrooms                   -17.799726
31          pipeline__total_rooms population                   -33.935309
33       pipeline__total_rooms median_income                   -11.407038
36       pipeline__total_bedrooms households                    32.750787
37    pipeline__total_bedrooms median_income                     7.849898
38                    pipeline__population^2                    27.458617
43                 pipeline__median_income^2                    24.379390
44      remainder__ocean_proximity_<1H OCEAN                    12.216885
45         remainder__ocean_proximity_INLAND                   -61.848075

Interpretabilidad en modelos no paramétricos

Los modelos de ensamble mejoran la precisión de las predicciones pero reducen su interpretabilidad por su aproximación de caja negra. En el ámbito del aprendizaje estadístico, se distinguen dos estrategias para abordar la interpretación de modelos predictivos:

• Enfoques globales: permiten interpretar la asociación entre predictores y respuesta al nivel del modelo, es decir, describen el comportamiento promedio del modelo para un conjunto de datos. Son enfoques globales:
  • Importancia de las variables por reducción en la función de coste (dependiente del modelo)
  • Importancia por permutación
  • Gráfico de dependencia parcial
• Enfoques locales: explican predicciones individuales, es decir, asignan a cada predicción individual la contribución de los diferentes predictores (si bien es posible construir explicaciones globales a partir de las distribuciones individuales).
  • Valores de Shapley

Conceptos

• Enfoque dependiente del modelo base: árboles
• Asocia la influencia de cada una de las variables del modelo predictivo a la reducción del criterio de generación de particiones (reducción de la impureza):
  • En clasificación: incertidumbre en un nodo (Gini o entropía).
  • En regresión: variabilidad de la respuesta dentro del nodo cuantificada por el Error Cuadrático (SSR).
    
• En concreto, la importancia de un predictor se calcula como la reducción total (normalizada) del criterio de división atribuible a esa característica:
  • Importancia del predictor \(j\) en árboles de clasificación: reducción en el índice de Gini a lo largo de los árboles por particiones sobre \(j\).
    
  • Importancia del predictor \(j\) en árboles de regresión: reducción en SSR a lo largo de los árboles por particiones sobre \(j\).

sklearn: los estimadores basados en árboles incluyen el atributo features_importances_

Limitaciones

• Al evaluarse sobre los datos de entrenamiento, puede sobrevalorar variables cuando el modelo sobreajusta asignando alta importancia a variables que no generalizan bien.

• Es una medida de la importancia sesgada: tiende a favorecer variables con alta cardinalidad (numéricas o con muchos valores posibles) frente a variables binarias o categóricas simples.

• Indica qué características son importantes en promedio, pero no explican su asociación con las predicciones del modelo: cómo afectan al resultado, es decir, dirección o forma de la relación.

Ejemplo: feature importance

from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_jobs=-1,max_features=0.5);
rf.fit(X_train,y_train);

# medidas de la importancia y ordenación
# en sklearn solo los árboles y ensembles de árboles (excepto bagging) tienen el artibuto feature_importances_
importances = rf.feature_importances_
indices = np.argsort(importances)

# representación gráfica
fig, ax = plt.subplots()
ax.barh(range(len(importances)), importances[indices])
ax.set_yticks(range(len(importances)))
_ = ax.set_yticklabels(np.array(X_train.columns)[indices])

Definición

Enfoque interpretativo independiente del modelo base. La importancia por permutación evalúa el aumento del error al permutar los valores de una variable en el conjunto de entrenamiento o de validación.

• Mide cuánto cambia el rendimiento del modelo cuando se permutan aleatoriamente (es decir, se “barajan”) los valores de una variable.

• La permutación puede hacerse:

  • dentro de una misma muestra, o
  • intercambiando valores entre muestras.

Si el modelo ha ajustado/aprendido una fuerte dependencia entre variable y respuesta, alterar esa variable cambia notablemente las predicciones y el rendimiento empeora.

Implementación

• Una variable \(k\) se considera importante si el error del modelo tras permutar \(k\), es decir \(e^{perm}_k\), es mucho mayor que el error original \(e\) \[ e^{perm}_k \gg e\]

• El valor de \(e^{perm}_k\) depende de la aleatorización de la permutación, y su variabilidad puede ser alta, sobre todo con muestras pequeñas. Para ello es usual:

  • repetir la permutación \(i\) veces,
  • promediar los errores resultantes, y
  • estimar la variabilidad muestral (incertidumbre) de la importancia.

Al romper la relación entre el predictor y la respuesta, se puede evaluar cuánto depende el modelo de ese predictor en particular. La importancia por permutación no refleja el valor predictivo intrínseco de una característica por sí sola, sino lo importante que es esa característica para un modelo en particular.

Importancia en muestra de entrenamiento y prueba

• La importancia por permutación pueden calcularse tanto en el conjunto de entrenamiento como en un conjunto de prueba o validación independiente.
  • Al calcular la importancia de permutación en el conjunto de entrenamiento identificamos qué variables contribuyen más al ajuste del modelo a los datos observados.
  • Usar un conjunto independiente permite identificar qué variables contribuyen más al poder de generalización del modelo analizado.
• La existencia de variables que resultan importantes en el conjunto de entrenamiento pero no en el de prueba puede ser indicativa de sobreajuste.

Ejemplo: importancia por permutación

import os
# La importancia de las características por permutación depende de la función de puntuación especificada mediante el argumento scoring (acepta múltiples funciones de puntuación). 
# Si no se indica ningún puntaje utiliza el criterio de evaluación por defecto del estimador. Ver detalles de la implementación en ayuda sklearn

from sklearn.inspection import permutation_importance
resultado = permutation_importance(rf, X_test, y_test,
                           n_repeats=30,
                           random_state=42,
                           scoring='r2',
                           n_jobs=os.cpu_count()-2)

# Los valores de la importancia por permutation nos indican la fraccion de caida del puntaje de referencia (R2) al permutar un predictor (destruir la info de esa variable)
#    Valores positivos y grandes: la variable es importante (su permutación empeora mucho el R2)
#    Valores cercanos a 0: poca o nula importancia (el modelo predice casi igual sin esa variable)
#    Valores negativos: la variable es perjudicial para el modelo (p.e. ruido). Al permutarla, el modelo mejora ligeramente
# IMPORTANTE: la interpretación es dependiente del scoring elegido

r2_test = rf.score(X_test,y_test)
print('R2 conjunto prueba:', r2_test)

df_imp = pd.DataFrame({
    "Feature": rf.feature_names_in_,
    "Importancia media": resultado.importances_mean,
    "Desviación": resultado.importances_std
}).sort_values("Importancia media", ascending=False)

df_imp

R2 conjunto prueba:0.812716732486059

sorted_importances_idx = resultado.importances_mean.argsort()
importancia = pd.DataFrame(
    resultado.importances[sorted_importances_idx].T,
    columns=X_train.columns[sorted_importances_idx],
)
ax = importancia.plot.box(vert=False, whis=10)
ax.set_title("Importancia de  permutación (conjunto test)")
ax.axvline(x=0, color="k", linestyle="--")
ax.set_xlabel("Disminución del score")
ax.figure.tight_layout()

Importancia por permutación: ventajas y limitaciones

• La importancia por permutación es una medida intuitiva y computacionalmente económica, ya que no requiere reentrenar el modelo, solo generar nuevas predicciones con los datos permutados.

• Sin embargo, puede ser inconsistente: si dos características contienen información similar, permutar una no refleja su verdadera importancia relativa. Solo permutar ambas o excluir una lo haría.

• Como medida global, la permutación por importancia solo indica qué variables son importantes, pero no cómo contribuyen al modelo (forma funcional ni dirección de la asociación).

Concepto

Los Valores de Shapley son un concepto prestado de la teoría de juegos cooperativos donde se usan para medir la distribución de un premio entre los distintos jugadores a partir de su contribución al mismo. Aplicado a un modelo, los valores de Shapley permiten atribuir a cada predictor la predicción de una observación específica.

• Se trata de una interpretación local (para una observación). En concreto la descomposión de Shapley \(\Phi^S(x_i;f)\) de un modelo \(f\) con \(p\) predictores para el la predicción evaluada en \(x_i\) es \[\Phi^S(x_i;f)=\phi_0^S+\sum_{k=1}^p\phi^s_k(x_i)=f(x_i)\]
  • \(\phi^S_0=E[f(X)]\) es el valor base o expectativa de la predicción del modelo (media de las predicciones del modelo)
  • \(\phi^S_k\) es la suma ponderada de la contribución marginal de la variable \(k\) para todas las posibles combinaciones de variables \(x\) al excluir la variable \(k\)-ésima (\(C(\mathbf{x})\, \backslash \, k\))
  • los valores de Shapley son aditivos: para la predicción en \(x_i\) suman \(f(x_i)\)
• Es posible construir una explicación global (a nivel de modelo) a partir de la distribución de los valores para las observaciones de entrenamiento o prueba. De la distribución se derivan tanto la importancia global de un predictor como su influencia en la dirección de cambio en las predicciones.

Valores de Shapley: estimación

• Sea una variable \(k\); la expresión \(\phi^s_k(x_i)\) evalúa cómo cambia la predicción en \(x_i\) al añadir o quitar una variable de distintos subconjuntos de variables. P.e., con tres variables \((x_1,x_2,x_3)\) que toman tres valores \(x_i=(A,B,C)\), la contribución de \(x_3=C\) sería \[\phi^S_C(f)=w_1[f(\{A,B,C\})-f(\{A,B\})]+w_2[f(\{A,C\})-f(\{A\})]+w_3[f(\{B,C\})-f(\{B\})]+w_4[f(\{C\})-f(\{\emptyset\})]\]

• El cálculo los valores de Shapley para todos los subconjuntos posibles de variables es computacionalmente muy costoso, ya que el número de subconjuntos crece exponencialmente con el número de variables.

• Para simplificar el cálculo, se utiliza muestreo Monte Carlo: se seleccionan aleatoriamente subconjuntos de variables en lugar de evaluarlos todos. En este caso se comparan las predicciones del modelo en subconjuntos aleatorios con y sin la variable \(k\) para luego promediar esas diferencias y estimar así la contribución de dicha variable.

• La desigualdad de Hoeffding proporciona una garantía estadística que acota la probabilidad de desviación entre los valores estimados y los reales: para un número suficientemente grande de muestras (coaliciones) \(K\), la probabilidad de que la estimación del valor de Shapley se desvíe del valor real más allá de un margen pequeño \(\epsilon\) es muy baja.

Implementación: el paquete shap

shap estima Valores de Shapley para explicar la salida de cualquier modelo de aprendizaje automático. Características:

• Una llamada típica consiste en instanciar un método o explainer y, a continuación, generar la explicación.
• Métodos soportados:
  • TreeExplainer (árboles/ensembles tipo XGBoost, LightGBM, RandomForest) con algoritmos rápidos y exactos/semiexactos.
  • KernelExplainer (válido para cualquier modelo) con aproximación por muestreo.
  • LinearExplainer (modelos lineales/GLM) eficiente.
  • DeepExplainer/GradientExplainer (redes profundas en TF/PyTorch) usa backprop/gradientes.
• Visualización:
  • Explicaciones individuales: force_plot, waterfall_plot.
  • Agregación global: summary_plot, bar_plot, beeswarm.
  • Visualización de efectos de variables sobre respuesta: scatter (forma y dirección del efecto).

Ejemplo: valores de Shapley

rf = RandomForestRegressor(n_estimators=100,
                            max_features='sqrt',
                            bootstrap=True,
                            n_jobs=6,
                            random_state=412)
rf.fit(X_train,y_train);

# Importamos libreria y creamos una instancia para explicar el modelo
import shap
explainer = shap.TreeExplainer(rf)
# Explainers específicos para distintos tipos de modelos (ver ayuda)
# A partir del explainer, explicamos las predicciones del modelo para el conjunto de datos de prueba
explicacion = explainer(X_test)

# La ejecución puede ser lenta. En concreto en un Mac M2 con 8gb de ram, la linea anterior
# se ejecuta en 12 min. Se puede optar por "explicar" un subconjunto aleatorio de 
# observaciones del conjunto de prueba. La linea a continuación se ejecuta en aprox 1:30 minutos
# explicacion = explainer(X_test.sample(500,random_state=42))

import matplotlib.pyplot as plt
fig = plt.figure()
shap.plots.waterfall(explicacion[1], show=False)
plt.gcf().set_size_inches(20/1.1,12/1.1)
plt.show()

fig = plt.figure()
shap.plots.waterfall(explicacion[2],show=False)
plt.gcf().set_size_inches(20/1.1,12/1.1)
plt.show()

fig = plt.figure()
shap.plots.beeswarm(explicacion, show=False)
plt.gcf().set_size_inches(20/1.1,12/1.1)
plt.show()

fig = plt.figure()
shap.plots.scatter(explicacion[:,'population'],color=explicacion,show=False)
plt.gcf().set_size_inches(20,12)
plt.show()

Inferencia: regresión de Shapley

A partir de los valores de Shapely podemos formular la regresión de Shapley \[y_i=\sum_{k=0}^p\phi_k^S(f,x_i) \beta_k^S+\hat{\epsilon_i}\qquad \hat{\epsilon_i}\sim N(0,\sigma^2_\epsilon)\] lo que nos permite estimar pseudo-coeficientes \(\beta_k^S\) y contrastar la hipótesis nula:

\[H_0^k(\Omega):\{ \beta_k^S\leq 0| \Omega\}\]

• La diferencia con el contraste habitual de hipótesis es su naturaleza local (\(\Omega\))
• Valores positivos indican variables significativas
• Los coeficientes \(\beta^S\) informan sobre la fuerza de la asociación entre la respuesta y el predictor \(k\).